/*=========================================================================

Program:   Medical Imaging & Interaction Toolkit
Language:  C++
Date:      $Date: 2010-07-31 05:33:05 +0800 (周六, 31 七月 2010) $
Version:   $Revision: 25026 $

Copyright (c) German Cancer Research Center, Division of Medical and
Biological Informatics. All rights reserved.
See MITKCopyright.txt or http://www.mitk.org/copyright.html for details.

This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE.  See the above copyright notices for more information.

=========================================================================*/

#if(_MSC_VER==1200)
#include <itkFixedCenterOfRotationAffineTransform.h>
#include <itkIndex.h>
#include <itkConceptChecking.h>
#endif
#include "WxRenderToContourCutter.h"

#include <mitkTimeHelper.h>
#include <mitkImageAccessByItk.h>
#include <mitkGeometry3D.h>
#include <vtkMatrix4x4.h>
#include <vtkCamera.h>
#include <vtkLinearTransform.h>

namespace mitk
{
    bool ipMITKSegmentationIsInsideContour( float *contour, int sizeContour, float x, float y )
    {

      int i, j;
      bool res = false;
      for ( i=0, j=sizeContour-1;
            i<sizeContour;
            j=i++)
      {
        if (
                (     ((contour[2*i+1]<=y) && (y<contour[2*j+1]))
                   || ((contour[2*j+1]<=y) && (y<contour[2*i+1]))  )
             &&
                (  x < (contour[2*j] - contour[2*i]) * (y - contour[2*i+1]) / (contour[2*j+1] - contour[2*i+1]) + contour[2*i])
           )

          res = !res;
      }
      return res;
    }

    template < typename TPixel, unsigned int VImageDimension, typename TOutputPixel >
    void CutImageWithOutputTypeSelect
      ( itk::Image<TPixel, VImageDimension>* inputItkImage, mitk::WxRenderToContourCutter* cutter, int /* boTimeStep */, TOutputPixel* /* dummy */)
    {
      typedef itk::Image<TPixel, VImageDimension> ItkInputImageType;
      typedef itk::Image<TOutputPixel, VImageDimension> ItkOutputImageType;
      typedef typename itk::ImageBase<VImageDimension>::RegionType ItkRegionType;
      typedef itk::ImageRegionIteratorWithIndex< ItkInputImageType > ItkInputImageIteratorType;
      typedef itk::ImageRegionIteratorWithIndex< ItkOutputImageType > ItkOutputImageIteratorType;

      if(cutter->m_OverlayContour.IsNull()||(cutter->m_Renderer==NULL))
        return;

      mitk::Contour::PointsContainerPointer points = cutter->m_OverlayContour->GetPoints();
      mitk::Contour::PointsContainerIterator pointsIt = points->Begin();

      float * contour=new float[points->Size()*2];
      int maxx=0,minx=32767,maxy=0,miny=32767;
      int i=0;
      Contour::BoundingBoxType::PointType point;
      while ( pointsIt != points->End() )
      {
        point = pointsIt.Value();
        contour[2*i]= point[0];
        contour[2*i+1]=point[1];
        if(point[0]>maxx) maxx=point[0];
        if(point[0]<minx) minx=point[0];
        if(point[1]>maxy) maxy=point[1];
        if(point[1]<miny) miny=point[1];

        pointsIt++;
        i++;
      }
      int contourSizeX=maxx-minx,contourSizeY=maxy-miny;
      bool *contourMask=new bool[contourSizeX*contourSizeY];

      for(int x=0;x<contourSizeX;x++)
      {
          for(int y=0;y<contourSizeY;y++)
          {
              if(ipMITKSegmentationIsInsideContour(contour,points->Size(),minx+x,miny+y))
                contourMask[y*contourSizeX+x]=true;
              else
                contourMask[y*contourSizeX+x]=false;
          }
      }
       delete[]  contour;

      if (inputItkImage == NULL)
      {
        return;
      }


      // first convert the index
      typename ItkRegionType::IndexType::IndexValueType tmpIndex[3];
      itk2vtk(cutter->m_InputRequestedRegion.GetIndex(), tmpIndex);
      typename ItkRegionType::IndexType index;
      index.SetIndex(tmpIndex);

      // then convert the size
      typename ItkRegionType::SizeType::SizeValueType tmpSize[3];
      itk2vtk(cutter->m_InputRequestedRegion.GetSize(), tmpSize);
      typename ItkRegionType::SizeType size;
      size.SetSize(tmpSize);

      //create the ITK-image-region out of index and size
      ItkRegionType inputRegionOfInterest(index, size);

      // PART 2: get access to the MITK output image via an ITK image
      typename mitk::ImageToItk<ItkOutputImageType>::Pointer outputimagetoitk = mitk::ImageToItk<ItkOutputImageType>::New();
      outputimagetoitk->SetInput(cutter->m_OutputTimeSelector->GetOutput());
      outputimagetoitk->Update();
      typename ItkOutputImageType::Pointer outputItkImage = outputimagetoitk->GetOutput();

      // PART 3: iterate over input and output using ITK iterators

      // create the iterators
      ItkInputImageIteratorType  inputIt( inputItkImage, inputRegionOfInterest );
      ItkOutputImageIteratorType outputIt( outputItkImage, outputItkImage->GetLargestPossibleRegion() );

      // Cut the boundingbox out of the image by iterating through
      // all pixels and checking if they are inside using IsInside()
      cutter->m_OutsidePixelCount = 0;
      cutter->m_InsidePixelCount = 0;
      mitk::Geometry3D* inputGeometry = cutter->GetInput()->GetGeometry();
      vtkMatrix4x4 *matrix = vtkMatrix4x4::New();
      vtkMatrix4x4 *IndexToWorldMatrix=vtkMatrix4x4::New();
      vtkMatrix4x4 *WorldToViewMatrix=vtkMatrix4x4::New();
      double view[4],Viewport[4];
      double dx,dy,dz;
      int renderWindowSize[2];
      renderWindowSize[0] = (cutter->GetRenderer()->GetRenderWindow()->GetSize())[0];
      renderWindowSize[1] = (cutter->GetRenderer()->GetRenderWindow()->GetSize())[1];
      cutter->GetRenderer()->GetViewport(Viewport);
      WorldToViewMatrix->DeepCopy(cutter->m_Renderer->GetActiveCamera()->
                       GetCompositeProjectionTransformMatrix(cutter->GetRenderer()->
                                                             GetTiledAspectRatio(),0,1));
      IndexToWorldMatrix->DeepCopy( inputGeometry->GetVtkTransform()->GetMatrix());
      vtkMatrix4x4::Multiply4x4(WorldToViewMatrix,IndexToWorldMatrix,matrix);
      IndexToWorldMatrix->Delete();
      WorldToViewMatrix->Delete();
      long time=0,time1=0;
      TOutputPixel outsideValue= (TOutputPixel) cutter->m_OutsideValue;
      //shall we use a fixed value for each inside pixel?

      clock_t tAll = clock();

      if (cutter->GetFillInside())
      {
        TOutputPixel insideValue  = (TOutputPixel) cutter->m_InsideValue;
        // yes, use a fixed value for each inside pixel (create a binary mask of the bounding object)

        for ( inputIt.GoToBegin(), outputIt.GoToBegin(); !inputIt.IsAtEnd(); ++inputIt, ++outputIt)
        {

//          clock_t t1 = clock();
          dx=inputIt.GetIndex()[0],dy=inputIt.GetIndex()[1],dz=inputIt.GetIndex()[2];
          view[0] = dx*matrix->Element[0][0] + dy*matrix->Element[0][1] +
            dz*matrix->Element[0][2] + matrix->Element[0][3];
          view[1] = dx*matrix->Element[1][0] + dy*matrix->Element[1][1] +
            dz*matrix->Element[1][2] + matrix->Element[1][3];
          view[2] = dx*matrix->Element[2][0] + dy*matrix->Element[2][1] +
            dz*matrix->Element[2][2] + matrix->Element[2][3];
          view[3] = dx*matrix->Element[3][0] + dy*matrix->Element[3][1] +
            dz*matrix->Element[3][2] + matrix->Element[3][3];
          if (view[3] != 0.0)
            {
            dx = view[0]/view[3];
            dy = view[1]/view[3];
            dz = view[2]/view[3];
            }
          dx = (dx + 1.0) * (renderWindowSize[0] *(Viewport[2]-Viewport[0])) / 2.0 + renderWindowSize[0]*Viewport[0];
          dy = (dy + 1.0) * (renderWindowSize[1] *(Viewport[3]-Viewport[1])) / 2.0 + renderWindowSize[1]*Viewport[1];

//          clock_t t2 = clock();
//          time+=(t2-t1);
          int xi=(int)dx-minx,yi=(int)dy-miny;
          if(xi>=0&&xi<contourSizeX&&yi>=0&&yi<contourSizeY&&contourMask[yi*contourSizeX+xi])
          {
            outputIt.Set(insideValue);
            ++cutter->m_InsidePixelCount;
          }
          else
          {
            outputIt.Set( (TOutputPixel) inputIt.Value());
            ++cutter->m_OutsidePixelCount;
          }
//           time1+=(clock()-t2);
        }
      }
      else
      {
        // no, use the pixel value of the original image (normal cutting)
        for ( inputIt.GoToBegin(), outputIt.GoToBegin(); !inputIt.IsAtEnd(); ++inputIt, ++outputIt)
        {

//          clock_t t1 = clock();
          dx=inputIt.GetIndex()[0],dy=inputIt.GetIndex()[1],dz=inputIt.GetIndex()[2];
          view[0] = dx*matrix->Element[0][0] + dy*matrix->Element[0][1] +
            dz*matrix->Element[0][2] + matrix->Element[0][3];
          view[1] = dx*matrix->Element[1][0] + dy*matrix->Element[1][1] +
            dz*matrix->Element[1][2] + matrix->Element[1][3];
          view[2] = dx*matrix->Element[2][0] + dy*matrix->Element[2][1] +
            dz*matrix->Element[2][2] + matrix->Element[2][3];
          view[3] = dx*matrix->Element[3][0] + dy*matrix->Element[3][1] +
            dz*matrix->Element[3][2] + matrix->Element[3][3];
          if (view[3] != 0.0)
            {
            dx = view[0]/view[3];
            dy = view[1]/view[3];
            dz = view[2]/view[3];
            }
          dx = (dx + 1.0) * (renderWindowSize[0] *(Viewport[2]-Viewport[0])) / 2.0 + renderWindowSize[0]*Viewport[0];
          dy = (dy + 1.0) * (renderWindowSize[1] *(Viewport[3]-Viewport[1])) / 2.0 + renderWindowSize[1]*Viewport[1];

//          clock_t t2 = clock();
//          time+=(t2-t1);
          int xi=(int)dx-minx,yi=(int)dy-miny;
          if(xi>=0&&xi<contourSizeX&&yi>=0&&yi<contourSizeY&&contourMask[yi*contourSizeX+xi])
          {
            outputIt.Set((TOutputPixel) inputIt.Value() );
            ++cutter->m_InsidePixelCount;
          }
          else
          {
            outputIt.Set( outsideValue);
            ++cutter->m_OutsidePixelCount;
          }
//            time1+=(clock()-t2);
        }
      }
      delete[] contourMask;
      matrix->Delete();

      long timeAll=clock()-tAll;

      MITK_INFO << "vtk_calc"<<time<<"other"<<time1<<"calc"<<timeAll;

    }

    template < typename TPixel, unsigned int VImageDimension >
    void CutImage( itk::Image< TPixel, VImageDimension >* inputItkImage, mitk::WxRenderToContourCutter* cutter, int boTimeStep )
    {
      TPixel* dummy = NULL;
      CutImageWithOutputTypeSelect<TPixel, VImageDimension, TPixel>(inputItkImage, cutter, boTimeStep, dummy);
    }

void WxRenderToContourCutter::SetOverlayContour( const mitk::Contour* overlayContour )
{
      m_OverlayContour = const_cast<mitk::Contour*>(overlayContour);
}

const mitk::Contour* WxRenderToContourCutter::GetOverlayContour() const
{
  return this->m_OverlayContour.GetPointer();
}

WxRenderToContourCutter::WxRenderToContourCutter()
  : m_OverlayContour(NULL), m_InsideValue(0), m_OutsideValue(0),m_Renderer(NULL),
    m_FillInside(true), m_OutsidePixelCount(0), m_InsidePixelCount(0)
{
  m_InputTimeSelector  = mitk::ImageTimeSelector::New();
  m_OutputTimeSelector = mitk::ImageTimeSelector::New();
}

WxRenderToContourCutter::~WxRenderToContourCutter()
{
}

const std::type_info& WxRenderToContourCutter::GetOutputPixelType()
{
  return *this->GetInput()->GetPixelType().GetTypeId();
}

void WxRenderToContourCutter::GenerateInputRequestedRegion()
{
  mitk::Image* output = this->GetOutput();
  if((output->IsInitialized()==false) )
    return;

  GenerateTimeInInputRegion(output, const_cast< mitk::Image * > ( this->GetInput() ));
}

void WxRenderToContourCutter::GenerateOutputInformation()
{
  mitk::Image::Pointer output = this->GetOutput();
  if ((output->IsInitialized()) && (output->GetPipelineMTime() <= m_TimeOfHeaderInitialization.GetMTime()))
    return;

  mitk::Image::Pointer input = const_cast< mitk::Image * > ( this->GetInput() );

  itkDebugMacro(<<"GenerateOutputInformation()");

  if(input.IsNull())
    return;  


  mitk::Geometry3D* inputImageGeometry = input->GetSlicedGeometry();
  m_InputRequestedRegion = input->GetLargestPossibleRegion();
  input->SetRequestedRegion(&m_InputRequestedRegion);


  unsigned int dimension = input->GetDimension();
  unsigned int *dimensions = new unsigned int [dimension];
  itk2vtk(m_InputRequestedRegion.GetSize(), dimensions);
  if(dimension>3)
    memcpy(dimensions+3, input->GetDimensions()+3, (dimension-3)*sizeof(unsigned int));
  output->Initialize(mitk::PixelType(GetOutputPixelType()), dimension, dimensions);
  delete [] dimensions;

  // now we have everything to initialize the transform of the output
  mitk::SlicedGeometry3D* slicedGeometry = output->GetSlicedGeometry();

  // set the transform: use the transform of the input; 
  // the origin will be replaced afterwards
  AffineTransform3D::Pointer indexToWorldTransform = AffineTransform3D::New();
  indexToWorldTransform->SetParameters(input->GetSlicedGeometry()->GetIndexToWorldTransform()->GetParameters());
  slicedGeometry->SetIndexToWorldTransform(indexToWorldTransform);
  
  // Position the output Image to match the corresponding region of the input image
  const mitk::SlicedData::IndexType& start = m_InputRequestedRegion.GetIndex();
  mitk::Point3D origin; vtk2itk(start, origin);
  inputImageGeometry->IndexToWorld(origin, origin);
  slicedGeometry->SetOrigin(origin);

  mitk::TimeSlicedGeometry* timeSlicedGeometry = output->GetTimeSlicedGeometry();
  timeSlicedGeometry->InitializeEvenlyTimed(slicedGeometry, output->GetDimension(3));
  timeSlicedGeometry->CopyTimes(input->GetTimeSlicedGeometry());

  m_TimeOfHeaderInitialization.Modified();
}

void WxRenderToContourCutter::ComputeData(mitk::Image* input3D, int boTimeStep)
{
   AccessFixedDimensionByItk_2(input3D, CutImage, 3, this, boTimeStep);
}

void WxRenderToContourCutter::GenerateData()
{
  mitk::Image::ConstPointer input = this->GetInput();
  mitk::Image::Pointer output = this->GetOutput();

  if(input.IsNull())
    return;  

  if((output->IsInitialized()==false) || (m_OverlayContour.IsNull()))
    return;

  m_InputTimeSelector->SetInput(input);
  m_OutputTimeSelector->SetInput(this->GetOutput());

  mitk::Image::RegionType outputRegion = output->GetRequestedRegion();
  const mitk::TimeSlicedGeometry *outputTimeGeometry = output->GetTimeSlicedGeometry();
  const mitk::TimeSlicedGeometry *inputTimeGeometry = input->GetTimeSlicedGeometry();

  ScalarType timeInMS;

  int timestep=0;
  int tstart=outputRegion.GetIndex(3);
  int tmax=tstart+outputRegion.GetSize(3);

  int t;
  for(t=tstart;t<tmax;++t)
  {
    timeInMS = outputTimeGeometry->TimeStepToMS( t );

    timestep = inputTimeGeometry->MSToTimeStep( timeInMS );

    m_InputTimeSelector->SetTimeNr(timestep);
    m_InputTimeSelector->UpdateLargestPossibleRegion();
    m_OutputTimeSelector->SetTimeNr(t);
    m_OutputTimeSelector->UpdateLargestPossibleRegion();
    clock_t t1 = clock();
        ComputeData(m_InputTimeSelector->GetOutput(), timestep);
    clock_t t2 = clock();
    long time=t2-t1;
     MITK_INFO << "all"<<time;
  }

  m_InputTimeSelector->SetInput(NULL);
  m_OutputTimeSelector->SetInput(NULL);

  m_TimeOfHeaderInitialization.Modified();
}

} // of namespace mitk
